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Abstract The Fast Marching Method is a very popular algorithm to com¬ 
pute times-of-arrival maps (distances map measured in time units). Since their 
proposal in 1995, it has been applied to many different applications such as 
robotics, medical computer vision, fluid simulation, etc. Many alternatives 
have been proposed with two main objectives: to reduce its computational 
time and to improve its accuracy. In this paper, we collect the main approaches 
which improve the computational time of the standard Fast Marching Method, 
focusing on single-threaded methods and isotropic environments. 9 different 
methods are studied under a common mathematical framework and experi¬ 
mentally in representative environments: Fast Marching Method with binary 
heap. Fast Marching Method with Fibonacci Heap, Simplified Fast March¬ 
ing Method, Untidy Fast Marching Method, Fast Iterative Method, Group 
Marching Method, Fast Sweeping Method, Lock Sweeping Method and Dou¬ 
ble Dynamic Queue Method. 

Keywords Fast Marching Method • Fast Sweeping Method • Fast Iterative 
Method • Eikonal Equation 


1 Introduction 

The Fast Marching Method and its derivatives (Fast Methods in short) have 
been extensively applied since they were firstly proposed in 1995 [35] as a 
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solution to isotropic control problems using first-order semi-Langragian dis¬ 
cretizations on Cartesians grids. Their main field of application are robotics 
[15,11,23] and computer vision [12], mainly medical image segmentation [4, 
2]. However, it has proven to be useful in many other applications such as 
tomography [26] or seismology [39]. 

The first approach was proposed by Tsitsiklis [35] but the most popular 
solution was given few months later by Sethian [29] using first-order upwind- 
finite differences in the context of isotropic front propagation. Differences and 
similarities between both works can be found in [33]. 

Fast Methods were originally proposed to simulate a wavefront propaga¬ 
tion through a regular discretization of the space. However, many different 
approaches have been proposed, extending these methods to other discretiza¬ 
tions and formulations. For a more detailed history of Fast Marching methods, 
we refer the interested readers to [7]. Despite the vast amount of Fast Methods 
literature, there is a lack of in-depth comparison and benchmarking. 

In this work 9 sequential (mono-thread), isotropic, grid-based Fast Meth¬ 
ods are detailed in the following sections: Fast Marching Method (FMM), 
Fibonacci-Heap FMM (FMMFib), Simplified FMM (SFMM), Untidy FMM 
(UFMM), Group Marching Method (GMM), Fast Iterative Method (FIM), 
Fast Sweeping Method (FSM), Lock Sweeping Method (LSM) and Double 
Dynamic Queue Method (DDQM). All these algorithms provide exactly the 
same solution except for UFMM and FIM, which have bounded errors. How¬ 
ever, the question of which one is the best for which applications is still open. 
For example, [16] compares only FMM and FSM in spite of the fact that 
GMM and UFMM were already published. This survey [20] mentions most 
of the algorithms but only compares FMM and SFMM. A more recent work 
compares FMM, FSM and FIM in 2D [6]. However, FIM is paralellized and 
implemented in CUDA providing a biased comparison. Fig. 1 schematically 
shows the comparisons among algorithms carried out in the literature. Meth¬ 
ods such as UFMM have been barely compared to their counterparts while 
others as FMM and FSM are compared in many works despite the fact that it 
is well known when each perform better: FSM is faster in simple environments 
with constant velocity. Also, results from one work cannot be directly extrap¬ 
olated to other works since the performance of these methods highly depend 
on their implementation. 

Statement of contributions: Three main contributions are included in 
this paper: 1) Based on previous works, a common formulation and notation is 
given for all the algorithms. This way, it is possible to easily understand their 
working principles and mathematical formulation. More detailed formulations 
are available in the literature but we decided to focus from a practical per¬ 
spective. 2) A recent survey of the work on designing sequential Fast Methods 
during the last years. We are not taking into account parallel or high-accuracy 
approaches as these fields are dense enough to fill another survey. 3) Extensive 
and systemic comparison among the mentioned methods with experiments de¬ 
signed taking into account their applications and results previously reported. 
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Fig. 1 Comparisons among algorithms. Colors refer to different works: orange [3], gray [16], 
yellow [19], green [6], black [37], and blue [20]. 


It is important to note that we have not taken into account the three meth- 
ods proposed in [7] (Fast Marching-Sweeping Method, Heap Cell Method and 
Fast Heap Cell Method) for 3 different reasons: 1) their error is not mathemat¬ 
ically bounded, 2) the error highly depends on the environment, and 3) they 
assume that the velocities are approximately the same in regions of arbitrary 
size, which is a strong assumption. 

This paper is organized as follows: next section introduces the common no¬ 
tation and formulation to be used throughout the document. Following, Fast 
Methods are detailed, classified attending to their algorithm family. Thus, Sec¬ 
tion 3 includes Fast Marching-like algorithms. Section 4 Fast Sweeping-based, 
and Section 5 contains other Fast Methods. The benchmark and its results 
are included in section Section 6, followed by a discussion in section Section 7. 
Finally, Section 8 outlines the conclusions and proposed future works. 


2 Problem Formulation 

Fast Methods are built to solve the nonlinear boundary value problem^. That 
is, given a domain Q and a function F : f2 ^ M_|_ which represents the local 
speed of the motion, drive a system from a starting set ffg C 12 to a goal set 
Xg C dl2 through the fastest possible path. The Eikonal equation computes 
the minimum time-of-arrival function T(x) as follows: 

|Vr(x)|F(x) = l,on 12 C . . 

T(x) = 0, X in dig ^ 

Once solved, T(x) represents a distances (times-of-arrival) field containing 
the time it takes to go from any point x to the closest point in d4 following 
the velocities on F(x). 

We assume, without loss of generality, that the domain is a unit hypercube 
of N dimensions: 12 = [0, l]'^. The domain is represented with a rectangular 
Cartesian grid X C containing the discretizations of the functions F(x) 
and T(x), F and T respectively. We refer to grid points Xjj = y^), Xjj S X 


^ This problem formulation closely follows [30] 
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as the point x = (cc, y) in the space corresponding to a cell (i, j) of the grid (for 
the 2D case). For notation simplicity, we will denote T-^ = T(xij) ss T(x), Tij € 
T, that is, Tjj represents an approximation to the real value of the function 
T(x). Analogously, Fy = F(xij) « F(x),Fij G F. For a general grid of N- 
dimensions, we will refer to cells by their index (or key) z as Xi, since a flat 
representation is more efficient for such datastructure. We will denote the set 
of Von-Neumann (4-connectivity in 2D) neighbors of grid point xy as A/'(xij). 


2.1 n-Dimensional Discrete Eikonal Equation 


In this section the most common first-order discretization of the Eikonal equa¬ 
tion is detailed. There exist many other first-order and higher-order approaches 
on grids, meshes and manifolds [22,32,1,24]. 

We are deriving the discrete Eikonal equation in 2D for better understand¬ 
ing. The most common first-order dicretization is given in [25], which uses 
an upwind-difference scheme to approximate partial derivatives of T(x) 
represents the one-sided partial difference operator in direction ±x): 


T,(x) 


Ty{^) 


pidzXrri _ Ti±ij ,Tij 

^ “ ±Zix 

r^^V'T' _ 


( 2 ) 


J max(DiT^T, 0)^ -p min(Dit^T, 0)^-^ 1 _ 1 
i max(Zljy^T, 0)^-I-min(D(|'^T, 0)^ J F)? 

Simpler but less accurate solution to Eq. (3) is proposed in [31]: 


r max(AT-T, - T, 0)^+ 1 _ J_ 

1 max(DiT*^T,-D,fT,0)2 j ^ > 

and Z\a; and Ay are the grid spacing in the x and y directions. Substituting 
Eq. (2) in Eq. (4) and letting 


T = r,.,- 

Tx = min(ri_i j, Tj+i j) (5) 

Ty = min(Tij_i,Tij+i) 

we can rewrite the Eikonal Equation, for a discrete 2D space as: 


max 


T-Tx 

/ix 


max 



1 


ij 


( 6 ) 


Since we are assuming that the speed of the front is positive (F > 0), 
T must be greater than Tx and Ty whenever the front wave has not already 
passed over the coordinates (z,j). Therefore, it is safe to simplify Eq. (6) to: 
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( t-t^ Y / r-Ty V ^ 

\ Ax ) \ Ay ) Fj? 


( 7 ) 


Eq. (7) is a regular quadratic equation of the form aT^ + 6T + c = 0, where: 


a = 
b = 


c = 


Z\2+Z\2 

-2{AlT^ + AlTy) 

Ain + AlT^-^ 


( 8 ) 


In order to simplify the notation for the n-dimensional case, let us assume 
that the grid is composed by cubic cells, that is, A^ = Ay = A^ = ■■■ = 
h. Let us denote Tj as the generalization of Tx or Ty for dimension d, up 
to N dimensions. We also denote as F the propagation velocity for point 
with coordinates (i,j, k,...). Therefore, the discretization of the Eikonal is a 
quadratic equation with parameters: 


a = N 

N 

b=-2J2Td 

d=l 


(9) 


2.2 Solving the nD discrete Eikonal equation 

Wavefront propagation follows causality. That is, in order to reach a point 
with higher time of arrival, it should firstly travel through neighbors of such 
point with smaller values. The opposite would imply a jump in time continuity 
and therefore the solutions would be erroneous. 

The proposed Eikonal solution (quadratic equation with parameters of 
Eq. (9)) does not guarantee the causality of the resulting distance map, as 
F and h can have arbitrary values. Therefore, before accepting a solution as 
valid its causality has to be checked. For instance, in 2D the Eikonal is solved 
as: 


Tx + Tj/ 1 2h^ , ,2 


( 10 ) 


called the two-sided update, as both parents Tx and Ty are taken into account. 
The solution is only accepted if T > max(Tx,ry). The upwind condition [7] 
shows that: 


T > max (Tx, Ty) 




If this condition fails, the one-sided update is applied instead: 


( 11 ) 
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r = min(T,,ry) + ^ (12) 

This is a top-down approach: the parents are iteratively discarded un¬ 
til a causal solution is found. To generalize Eq. (11) is complex. Therefore, 
we choose to use a bottom-up approach: Eq. (12) is solved and parents are 
iteratively included until the time of the next parent is higher than the cur¬ 
rent solution: Tk > T. The procedure is detailed in Algorithms 1 and 2. The 
MinTDim() function returns the minimum time of the neighbors in a given 
dimension (left and right for dim = 1, bottom and top for dim = 2, etc.). Our 
experiments found this approach more robust for 3 or more dimensions with 
negligible impact on the computational performance. 


Algorithm 1 Solve Eikonal Equation 
1: procedure SOLVEElKONAL(xi, T, 

2: N 

3: for dim = I : N do 

4: miriT <— MlNTDlM(dim) 

5: if minj- ^ oo and minx < T then 

6: Tvaiues -pushlminr) 

7: else 

8: a <— a — 1 

9: if a = 0 then t> FSM can cause this situation. 

10: return oo 

11: "lvalues SORT(Tvalues) 

12: for dim = 1 : a do 

13: Ti <- SOLVENDlMS(xi,(iim,Tvaiues,-F) 

14: if dim = a or Ti < 7^aiues,dim+l then 

15: break 

16: return Tj 


Algorithm 2 Solve Eikonal for n dimensions 

1: procedure SOLVENDlMS(xi, dim, Tvaiues, 1^) 

2: if dim = 1 then 

3: return Tvaiues.l + jk 

dim 

4. SUTTiT ^ y ^ ^^alues,! 

i=l 

dim 

5: sumT^ ^ E Thames,i 

i = l 

6: a ^ dim 

7: b < - 2sumT 

8: c ■<— sumT^ — ^ 

9: q — Aac 

10: if < 0 then > Non-causal solution 

11: return oo 

12: else 

13: return ->>+‘^^< 1 ) 

2a 
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3 Fast Marching Methods 

The Fast Marching Method (FMM) [29] is the most common Eikonal solver. 
It can be classified as a label-setting, Dijkstra-like algorithm [10]. It uses a 
first-order upwind-finite difference scheme to simulate an isotropic front prop¬ 
agation. The main difference with Dijsktra’s algorithm is the operation carried 
out on every node. Dijkstra’s algorithm is designed to work on graphs. There¬ 
fore, the value for every node Xi only depends on one parent xj, following the 
Bellman’s optimality principle [5]: 


Ti = min (cijTTj) 

XieV(xi) 


(13) 


In other words, a node Xi is connected to the parent xj in its neighborhood 
A/'(xi) which minimizes (or maximizes) the function value (in this case Ti) 
composed by the value of Tj plus the addition of the cost of traveling from xj 
to Xi, represented as c^. 

The FMM follows Bellman’s optimality principle but the value for every 
node is computed following first-order upwind discretization of the Eikonal 
equation detailed in Section 2. This discretization takes into account the spatial 
representation (i.e. a rectangular grid) and the value of all the causal upwind 
neighbors. Thus, the times-of-arrival field computed by FMM is more accurate 
than Dijkstra’s. 

The algorithm divides the cells in three different sets: I) Frozen: those 
cells which value is computed and cannot change, 2) Unknown: cells with no 
value assigned, to be evaluated, and 3) Narrow band (or just Narrow): frontier 
between Frozen and Unknown containing those cells with a value assigned 
that can be improved. These sets are mutually exclusive, that is, a cell cannot 
belong to more than one of them at the same time. The implementation of 
the Narrow set is a critical aspect of FMM. A more detailed discussion will be 
carried out in Section 3.1. 

The procedure is detailed in Algorithm 3. Initially, all points^ in the grid 
belong to the Unknown set with infinite arrival time. The initial points (wave 
sources) are assigned a value 0 and introduced in Frozen (lines 2-7). Then, the 
main FMM loop starts by choosing the element with minimum arrival time 
from Narrow (line 9). All its non-Frozen neighbors are evaluated: for each 
of them the Eikonal is solved and the new arrival time value is kept if it is 
improved. In case the cell is in Unknown, it is transferred to Narrow (lines 10- 
16). Finally, the previously chosen point from Narrow is transferred to Frozen 
(lines 17 and 18) and a new iteration starts until the Narrow set is empty. The 
arrival times map T is returned as the result of the procedure. 


From now on, we will indistinctly use point, cell or node to refer to each element of the 
grid. 
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Algorithm 3 Fast Marching Method 


1 

procedure FMM(A', T, T, A's) 
Initialization: 


2 

Unknown ^ Af, Narrow ^ 0, Frozen -f- 0 


3 

Ti ■<— oo Vxi G X 


4 

for Xi G Xs do 


5 

Ti ^ 0 


6 

Unknown Unknown\{xi} 


7 

Narrow •<— Narrow U {xi} 



Propagation: 


8 

while Narrow 7^ 0 do 


9 

Xmin t- argminxieMarrou {71} 

> Narrow top operation. 

10 

for x; g (Af(xmin) n A’\Frozen) do 

> For all neighbors not in Frozen. 

11 

Tj <- SOLVEElKONAL(xi,T, J') 


12 

if Ti < Ti then 


13 

Ti Ti 

> Narrow increase operation if xi G Narrow. 

14 

if Xj G Unknown then 

> Narrow push operation. 

15 

Narrow ^ Narrow U {xi} 


16 

Unknown Unknown\{xi} 


17 

Narrow Narrow\{xniin} 

> Narrow pop operation: add to Frozen. 

18 

Frozen Frozen U {xmin} 


19 

return T 



3.1 Binary and Fibonacci Heaps 

FMM requires the implementation of the Narrow set to have four different 
operations: 1) Push: to insert a new element to the set, 2) Increase: to reorder 
an element already existing in the set which value has been improved, 3) Top: 
retrieve the element with minimum value, and 4) Pop: remove the element 
with minimum value. As stated before, this is the most critical aspect of the 
FMM implementation. The most efficient way to implement Narrow is by using 
a min-heap data structure. A heap is an ordered tree in which every parent is 
ordered with respect to its children. In a min-heap, the minimum value is at 
the root of the tree and the children have higher values. This is satisfied for 
any parent node of the tree. 

Among all the existing heaps, FMM is usually implemented with a binary 
heap [17]. However, the Fibonacci Heap [13] has a better amortized time for 
Increase and Push operations. However, it has additional computational over¬ 
head with respect other heaps. For relatively small grids, where the narrow 
band is composed by few elements and the performance is still far from its 
asymptotic behavior, the binary heap performs better. Table 1 summarizes 
the time complexities for these heaps^ (the priority queue will be detailed in 
Section 3.2). Note that n for the methods complexity is the number of cells in 
the map, as the worst case is to have all the cells in the heap. 

Each cell is pushed and popped at most once in the heap. For each loop, the 
top of Narrow is accessed (0(1)), the Eikonal is solved for at most 2^ neighbors 

^ http://bigocheatsheet.com/ 
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Table 1 Summary of amortized time complexities for common heaps used in EMM (n is 
the number of elements in the heap). 



Push 

Increase 

Top 

Pop 

Fibonacci 

0(1) 

0(1) 

0(1) 

O(logn) 

Binary 

O(logn) 

O(logn) 

0(1) 

O(logn) 

Priority Queue 

O(logn) 

- 

0(1) 

O(logn) 


(0(1) for a given N), these cells are pushed or increased (O(logn) in the worst 
case), and finally the top cell is popped (O(logn)). Therefore each loop is at 
most O(logn). Since this loop is executed at most n times, the total FMM 
complexity is O(nlogn), where n represents the total number of cells of the 
grid in the worst case scenario. 


3.2 Simplified Fast Marching Method 

The Simplified Fast Marching Method (SFMM) [20] is a relatively unknown 
variation of the standard FMM but with an impressive performance. SFMM 
is a reduced version of FMM, where Narrow, implemented as a simple priority 
queue, can contain different instances of the same cell with different values. 
Additionally, it can happen that the same cell belongs to Narrow and Frozen 
at the same time. The simplification occurs since no Increase operation is 
required. Every time a cell has an updated value, it is pushed to the priority 
queue. Once it is popped and inserted in Frozen, the remaining instances in 
the queue are simply ignored. 

The advantage is that all the increase operations are substituted by push 
operations. Although both have the same computational complexity, the con¬ 
stant for push is much lower (increase requires removal and Push operations). 
Note that the computational complexity is maintained, O{n\ogn). 


3.3 Untidy Fast Marching Method 

The Untidy Fast Marching Method (UFMM) [37,27] follows exactly the same 
procedure as FMM. However, a special heap structure is used which reduces 
the overall computational complexity to 0{n): the untidy priority queue. 

This untidy priority queue is closer to a look-up table than to a tree. It 
assumes that the iF values are bounded, hence the T values are also bounded. 
The untidy queue is a circular array which divides the maximum range of T 
into a set of k consecutive buckets. Each bucket contains an unordered list of 
cells with similar Ti value. The threshold values for each bucket evolve with 
the algorithm, trying to maintain an uniform distribution of the elements in 
Narrow among the buckets. 

Since the index of the corresponding bucket can be analytically computed, 
Push is 0{1) as well as Top and Pop. The Increase operation is, in average, 
0{1) as long as ^buckets <0{n). Therefore, the total UFMM complexity is 
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Algorithm 4 Simplified Fast Marching Method 

1 

procedure SFMM(A', T, A's) 

Initialization as FMM (in Algorithm 3) 



Propagation: 


2 

while Narrow 7 ^ 0 do 


3 

Xmin <- argminxi6Karro» {71} 

> Narrow top operation. 

4 

if Xmin ^ Frozen then 


5 

Narrow •<— Narrow\{xinin} 


6 

else 


7 

for Xi g (Aflxinin) n A’\Frozen) do 

> For all neighbors not in Frozen. 

8 

Ti <- SOLVEElKONAL(xi,T, J") 


9 

if T[ < Tj then 

> Update arrival time. 

10 

71 <- Ti 


11 

Narrow Narrow U {xj} 

> Narrow push operation. 

12 

if Xi G Unknown then 


13 

Unknown Unknown\{xi} 


14 

Narrow Narrow\{xniin} 

> Narrow pop operation. 

15 

Frozen Frozen U {xmin} 


16 

return T 



0{n). However, since elements within a bucket are not sorted (FIFO strategy 
is applied in each bucket), errors are being introduced in the final result. In 
fact, it is shown that the accumulated additional error is bounded by 0{h), 
which is the same order of magnitude as in the original FMM. 


dt 


2dt 


cl c3 


c2 


(k-1)dt k dt (max. range) 

c4 


dt 


2dt 


3dt 


kdt (k+1)dt 


c2 c8 c7 


c6 c5 


(k-1)dt 


k-dt 


(k+1)dt 


c4 c36 


c42 c21 c84 


2(k-1)dt i 2kdt 
c92 c22 


Fig. 2 Untidy priority queue representation. Top: first iteration, the four neighbors of 
the initial point are pushed. Middle: the first bucket became empty, so the circular array 
advances one position. Cell c2 is first evaluated because it was the first pushed in the bucket. 
Bottom: after a few iterations, a complete loop on the queue is about to be completed. 
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4 Fast Sweeping Methods 

The Fast Sweeping Method (FSM) [34,40] is an iterative algorithm which 
computes the times-of-arrival map by successively sweeping (traversal) the 
whole grid following a specific order. FSM performs Gauss-Seidel iterations 
in alternating directions. These directions are chosen so that all the possible 
characteristic curves of the solution to the Eikonal are divided into the possible 
quadrants (or octants in 3D) of the environment. For instance, a bi-dimensional 
grid has 4 possible Gauss-Seidel iterations (the combinations of traversing 
X and y dimensions forwards and backwards): are North-East, North-West, 
South-East and South-West, as shown in Fig. 3. 



Fig. 3 FSM sweep directions in 2D represented with arrows. The darkest cell is the initial 
point and the shaded cells are those analyzed by the current sweep (time improved or 
maintained). 


The FSM is a simple algorithm: it performs sweeps until no value is im¬ 
proved. In each sweep, the Eikonal equation is solved for every cell. However, 
to generalize this algorithm to iV-dimensions is complex. Up to our knowledge, 
there are only 2D and 3D implementations. However, in Algorithm 5 we intro¬ 
duce an V-dimensional version. We will denote the sweeping directions as a 
binary array SweepDirs with elements 1 or —1, with 1 (—1) meaning forwards 
(backwards) traversal in that dimension. This array is initialized to 1 (North- 
East in the 2D case or North-East-Top in 3D) and the grid is initialized as in 
FMM (lines 2-5). The main loop updates SweepDirs and a sweep is performed 
in the new direction (lines 8-9). 

The GetSweepDirs() procedure (see Algorithm 6) is in charge of gener¬ 
ating the appropriate Gauss-Seidel iteration directions. If a 3D SweepDirs = 
[1,1,1] vector is given, the following sequence will be generated: 


1- [-l,-l,-l] 

2 - [ 1 ,- 1 ,- 1 ] 

3 -[- 1 , 1 ,- 1 ] 


4 - [ 1 , 1 ,- 1 ] 

5- [-l,-l,l] 

6- [l,-l,l] 


7 - [- 1 , 1 , 1 ] 

8 - [ 1 , 1 , 1 ] 


(14) 


Note that this sequence creates a sweep pattern which is not exactly the 
same as detailed in the literature, but it is equally valid as the same directions 
are visited and the same number of sweeps are required to cover the whole 
grid. 
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Finally, the Sweep() procedure (see Algorithm 7) recursively generates the 
Gauss-Seidel iterations following the traversal directions specified by the cor¬ 
responding value of SweepDirs (Line 4). Once the most inner loop is reached, 
the corresponding cell is evaluated and its value updated if necessary (lines 
7-11). 


Algorithm 5 Fast Sweeping Method 

1: procedure FSM(A’,r,J',A’s) 

Initialization. 

2: SweepDirs [1,..., 1] > Initialize sweeping directions. 

3: Ti oo Vxi G A' 

4: for G A's do 

5: Ti ^ 0 

Propagation: 

6 : stop False 

7: while stop 7 ^ True do 

8 : SweepDirs •<— GEtSweepDirs(A:’, SweepDirs) 

9: stop Sweep(A:', T, SweepDirs, A^) 

10 : return T 


Algorithm 6 Sweep directions algorithm 

1: procedure GEtSweepDirs(A’, SweepDirs) 

2: for i = 1 : N do 

3: SweepDirSj SweepDirs^ + 2 

4: if SweepDirSj < 1 then 

5: break > Finish For loop. 

6 : else 

7: SweepDirSj < -1 

8 : return SweepDirs 


Algorithm 7 Recursive sweeping algorithm 

1: procedure Sweep(A:’, T, 7^, SweepDirs, n) 

2: stop •<— True 

3: if n > 1 then 

4: for i G A’n following SweepDirs^^ do 

5: stop Sweep(A:’, T, 7^, SweepDirs, n — 1) 

6 : else 

7: for i G A’l following SweepDirS;!^ do 

8: Ti ■<— SOLVEEiKONAL(xi, T, ^) l> Xi is the corresponding cell. 

9: if T[ < Tj then 

10; Ti ^ Ti 

11: stop False 

12 : return stop 
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The FSM carries out as many grid traversals as necessary until the value 
Ti for every cell has converged. Since no ordering is implied, the evaluation of 
each cell is 0(1). As there are n cells, the total computational complexity of 
FSM is 0{n). However, note that the constants highly depend on the velocity 
function A(x). In the case of an empty environment with constant F{x), only 
2^ sweeps will be required as the characteristic directions are straight lines. 
However, for environment with obstacles or complex velocities functions, where 
the characteristics directions change frequently, the number of sweeps required 
can be much higher and therefore FSM will take longer to return a solution. 
Note as well that the T returned by FSM is exactly the same as all the FMM- 
like algorithms (except UFMM). 


4.1 Lock Sweeping Methods 

The Lock Sweeping Method (LSM) [3] is a natural improvement over FSM. 
The FSM might spend computation time recomputing T, even if none of the 
neighbors of x, has improved their value since the last sweep. LSM labels a 
cell as unlocked if any of its neighbors has changed and thus its value can be 
improved. Otherwise, the cell is labeled as locked and it will be skipped. 

The LSM procedure is detailed in Algorithm 8. It is basically the same 
as FSM but with the addition of tracking if a point is locked (Frozen) or un¬ 
locked (Narrow). Analogously, the LockSweep() (see Algorithm 9) procedure 
is similar to Sweep() with two differences: 1) if a point is not unlocked it is 
skipped (see Line 8), and 2) neighbors of Xi are unlocked if the new value Ti 
is better that their current value. 


Algorithm 8 Lock Sweeping Method 

1: procedure LSMIA”, T, T, Aij) 

Initialization. 

2: Frozen <—A”, Narrow <—0 

3: T; <— oo Vx; £ X 

4: SweepDirs <— [1,..., 1] > Initialize sweeping directions. 

5: for Xi S Xb do 

6: T; ^ 0 

7: for Xj £ Affxi) do l> Unlocking neighbors of starting cells. 

8: Frozen <—Frozen\{xj} 

9: Narrow <— Narrow U {xj } 

Propagation: 

10: stop <— False 

11: while stop ^ True do 

12: SweepDirs <— GETSwEEPDlRSfA”, SweepDirs) 

13: stop -f— LoCKSwEEPfA”, T, T, SweepDirs, TV) 

14: return T 


Note that the asymptotic computational complexity of FSM is kept, 0{n). 
The number of required sweeps is also maintained. However, in practice it 






14 


Javier V. Gomez et al. 


Algorithm 9 Recursive sweeping algorithm 
1: procedure LockSweep(A’, T, SweepDirs, n) 

2: stop i— True 

3: if n > 1 then 

4: for i G following SweepDirSj^ do 

5: stop LockSweep(A:’, T, SweepDirs, n — 1) 

6 : else 

7: for i G following SweepDirS;!^ do 

8 : if Xj G Narrow then 

9: Ti ■<— SOLVEEiKONAL(xi, T, 7^) t> Xi is the corresponding cell. 

10 ; if Ti < Ti^then 

11; Ti^Ti 

12; stop False 

13; for Xj G A/^(xi) do 

14; if Tj < Tj then > Add improvable neighbors to Narrow. 

15; Frozen ^ Frozen\{xj} 

16; Narrow ^ Narrow U {xj} 

17; Narrow Narrow\{xi} 0 Add Xj to Frozen. 

18; Frozen •<— Frozen U {xj} 

19; return stop 


turns out that most of the cells are locked during a sweep. Therefore, the 
computation time saved is important. 


5 Other Fast Methods 

5.1 Group Marching Method 

The Group Marching Method (GMM) [21] is an FMM-based Eikonal solver 
which solves a group of grid points in Narrow at once, instead of sorting them 
in a heap structure. 

Consider a front propagating. At a given time. Narrow will be composed 
by the set of cells belonging to the wavefront. GMM selects a group G out 
of Narrow composed by the global minimum and the local minima in Narrow. 
Then, every neighboring cell to G is evaluated and added to Narrow. These 
points in G have to be chosen carefully so that causality is not violated since 
GMM does not sort the Narrow set. For that, GMM selects those points fol¬ 
lowing Eq. (15): 


G = {xi G Narrow : Ti < min(TNarrow) + Sr} 


(15) 


where 




1 


(16) 


max(T') 

In the original GMM work [21], 6r = -However, we have chosen 

Eq. (16) as referred in [19]. Although this second formula is not mathematically 
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proven, the results for the original Sr are much worse than FMM in most of 
the cases, reaching an order of magnitude of difference. 

If the time difference between two adjacent cells is larger than Sr, their 
values will barely affect each other since the wavefront propagation direction 
is more perpendicular than parallel to the line segment formed by both cells. 
However, the downwind points (those to be evaluated in future iterations) can 
be affected by both adjacent cells. Therefore, points in G are evaluated twice 
to avoid instabilities. 

GMM is detailed in Algorithm 10. Its initialization is FMM-like. Note that 
Sr depends on the maximum velocity value in the grid. The main loop updates 
the threshold Tm every iteration. Firstly, it carries out a reverse traversal 
through the selected points, computing and updating their value (lines 17-21). 
Then, lines 22-31 perform a forward traversal with the same operations as the 
reverse traversal but updating the Narrow and Frozen sets in the same way 
as FMM. 

Note that GMM returns the same solution as FMM. GMM evaluates twice 
every node before inserting it in Frozen while FMM only evaluates it once. 
However, GMM does not require any sorting. Therefore, GMM is an 0(n) it¬ 
erative algorithm that converges in only 2 iterations (traversals). The value 
of Sr can be modified: higher Sr would require more iterations to converge. 
However, smaller Sr will require also 2 traversals but the group G will be com¬ 
posed by less cells. As GMM authors point out, GMM can be interpreted as 
an intermediary point between FMM (Sr = 0) and a purely iterative method 
[28] (Sr = oo). 

Additionally, a generalized A^-dimensional implementation is straightfor¬ 
ward. 


5.2 Dynamic Double Queue Method 

The Dynamic Double Queue Method (DDQM) [3] is inspired in the LSM but 
resembles to GMM. DDQM is conceptually simple. Narrow is divided into 
two non-sorted FIFO queues: one with cells to be evaluated sooner and the 
other one with cells to be evaluated later. Every iteration takes an element 
from the first queue and evaluates it. If the time is improved, the neighboring 
cells with higher time are unlocked and added to the first or second queue 
depending on the value of the cell updated. Once the first queue is empty, 
queues are swapped and the algorithm continues. The purpose is to achieve a 
pseudo-ordering of the cells, so that cells with lower value are evaluated first. 

Since queues are not sorted, it could require to solve many times the same 
cell until its value converges. DDQM dynamically computes the threshold value 
depending on the number of points inserted in each queue, trying to reach 
an equilibrium. The original paper includes an important analysis about this 
threshold update. Initially, the step value of the threshold is increased every 
iteration and is computed as: 
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Algorithm 10 Group Marching Method 


1 : procedure GMM(X,T, T, X^) 
Initialization: 


Unknown <— X, Narrow <— 0, Frozen 0 

T; <— OO Vx; S X 

^ max(J^) 

for Xi S A’s do 
Ti ^ 0 

Unknown <— Unknown\{xi} 

Frozen <— Frozen U{x;} 

for Xj £ A/'(xi) do O Adding neighbors of starting points to Narrow. 

Tj <- SOLVEElKONAL(xj, T, J') 
if Ti < Tm then 
Tm <- Ti 

Unknown <— Unknown\{xi} 

Narrow <— Narrow U {x;} 
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Propagation: 

while Narrow ^ 0 do 

Tm I’m + i5t- 

for Xi S (Narrow < Tm) REVERSE do 
for Xj S (A/’(xi) n A’\Frozen) do 
Ti <- SOLVEElKONAL(xj,T, T) 
if Ti < Ti then 
Ti <- Ti 

for Xi £ (Narrow < Tm) FORWARD do 
for Xj £ (A/’(xi) n A’\Frozen) do 
Ti •(- SOLVEElKONAL(xj,T, T) 
if Ti < Ti then 
Ti ^ Ti 

if Xi £ Unknown then 

Unknown <— Unknown\{xi} 
Narrow <— Narrow U {xi} 
Narrow <— Narrow\{xi} 

Frozen <— Frozen U {xi} 
return T 


> Reverse traversal. 


> Eorward tranversal. 


step 


1.5hn 


EPi 


(17) 


where n is the total number of cells in the grid. Originally, this step was 
proposed as step = j_ . However, step should have time units and this 

expression (probably an error due to the ambiguity of using speed F or 
slowness f = p)- Therefore, we propose as alternative Eq. (17). 

Every time the first queue is empty, UpdateStep() (see Algorithm 11) is 
called, with the value of the current step, ci, and Ctotai which are the number 
of cells inserted in the hrst queue and the total number of cells, correspond¬ 
ingly. step is modified so that the number of cells inserted in the hrst queue 
are between 65% and 75%. This is a conservative approach, since the closer 
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this percentage is to 50% the faster DDQM is. However, the penalization for 
percentages lower than 50% is much important than for higher percentages. 

Note that the step is increased by a factor 1.5 but decreased by a factor of 
2. This makes step to converge to a value instead of overshooting around the 
optimal value. Dividing by a larger number causes the first queue to become 
empty earlier. Thus, next iteration will finish faster and a better step value 
can be computed. 


Algorithm 11 DDQM Threshold Increase 
1: procedure UPDATESTEP(step, ci, Ctot^l) 

2: m ^ 0.65 

3: M -s- 0.75 

4: Perc <— 1 

5: if Cl > 0 then 

6 : Perc <— —21— 

^total 

7: if Perc < m then 

8 : step ■<— step * 1.5 

9: else if Perc > M then 

10 : step e- ^ 

11 : return step 


DDQM is detailed in Algorithm 12. As in LSM, points are locked (Frozen) 
or unlocked (Narrow). Initialization sets all points as frozen except the neigh¬ 
bors of the start points, which are added to the first queue (lines 2-14). While 
first queue is not empty, its front element is extracted and evaluated (lines 
16-18). If its value is improved, all its locked neighbors with higher value are 
unlocked and added to its corresponding queue. 

In the original work, three methods were proposed: 1) single-queue (SQ), 
and therefore simpler algorithm, 2) two-queue static (TQS), where the step 
is not updated, and 3) two-queue dynamic (which we call DDQM). SQ and 
TQS slightly improve DDQM in some experiments, but when DDQM improves 
SQ and TQS (for instance environments with noticeable speed changes) the 
difference can reach one order of magnitude. Therefore, we decided to include 
DDQM instead of SQ and TQS since it has shown a more adaptive behaviour. 
In any case, any of these methods return the same solution as FMM. 

In the worst case, the whole grid is contained in both queues and traversed 
many times during the propagation. However, since queue insertion and dele¬ 
tion are 0{1) operations, the overall complexity is 0(n). Note that SWAP() 
can be efficiently implemented in 0(1) as a circular binary index, or updating 
references (or pointers). There is not need for a real swap operation. 


5.3 Fast Iterative Method 

The Fast Iterative Method (FIM) [19] is based on the iterative method pro¬ 
posed by [28] but inspired in FMM. It also resembles to DDQM (concretely to 
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Algorithm 12 Double Dynamic Queue Method 

1: procedure DDQK'hA’, T, A’s) 

Initialization: 

2: Frozen <—A”, Narrow <—0 

3: Qi^0,S2^0 

4: Cl <— 0 > Counters 

5: Ctotal 0 

6: step = > n is the total number of cells. 

7: t/i <— step 

8: Ti <— oo Vxi £ X 

9: for x; £ X^ do 

10: Ti ^ 0 

11: for Xj £ A/^(xi) do 

12: Ql <— Qi U {xj} 

13: Unknown <— Unknown\{xi} 

14: Narrow <— Narrow U {xi} 
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Propagation: 

while Ql ^ 0 or Q 2 7 ^ 0 do 
while Ql 7 ^ 0 do 

Xi •«— Qi.fhont() > Extracts the front element. 

Ti <- SOLVEElKONAL(xi,T, J") 
if Ti < Ti then 
Ti t— Ti 

for Xj £ (A/'(xi) n Frozen) do 

if Ti < Tj then > Add improvable neighbors to corresponding queue. 
Frozen £- Frozen\{xj} 

Narrow •<— Narrow U {xj} 

Ctotal Ctotal + 1 

if Ti < th then 

Ql ^ Ql u {xj} 

Cl <- Cl + 1 
else 

Q2 ^ (Q2 U {xj}) 

Narrow •<— Narrow\{xi} 

Frozen Frozen U {xi} 
step ^ UPDATESTEP(step, Cl, Ctotal) 
swap(Qi, Q2) 

Cl i — 0 
Ctotal ■<— 0 
th ^ th step 

return T 


its single queue variant). It iteratively evaluates every point in Narrow until it 
converges. Once a node has converged its neighbors are inserted into Narrow 
and the process continues. Narrow is implemented as a non-sorted list. The 
algorithm requires a convergence parameter e: if Ti is improved less than e, it 
is considered converged. FIM has been also proposed for triangulated surfaces 
[14]. 

FIM is designed to be efficient for parallel computing, since all the ele¬ 
ments in Narrow can be evaluated simultaneously. However, we are focusing 
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on its sequential implementation in order to have a fair comparison with other 
methods. 

Algorithm 13 details FIM. Its initialization is the same as FMM. Then, 
for each element in Narrow, its value is updated (lines 11-12). If the value 
difference is less than e, the neighbors are evaluated and added to Narrow 
in case their value is improved (lines 13-19). Since Narrow is a list, the new 
elements should be inserted just before the point being currently evaluated, 
Xi. Finally, this point is removed from Narrow and labeled as Frozen (lines 20 
and 21). 


Algorithm 13 Fast Iterative Method 

1: procedure FIM(A', T, .F, A's, e) 

Initialization: 

2: Frozen ^ T, Narrow •<— 0 

3: Ti ■<— oo Vxi G X 

4: for Xi G Xs do 

5: Ti ^ 0 

6 : for Xj G (A/*{xi) H Unknown) do 

7: Frozen Frozen\{xi} 

8 : Narrow Narrow U {xi} 

Propagation: 

9: while Narrow 7 ^ 0 do 

10: for Xi G Narrow do 

11: Ti ^ Ti 

12: Ti SOLVEElKONAL(xi,T, T) 

13: if |Ti — Ti I < e then 

14: for Xj G (A/*(xi) H Frozen) do 

15: Tj SOLVEElKONAL(xj,T, T) 

16: if Tj < Tj then 

17: Tj ^ fj 

18: Frozen <—Frozen\{xj} 

19: Narrow <— Narrow U {xj} t> Insert in the list just before x; 

20: Narrow <— Narrow\{xi} 

21: Frozen <—Frozen U {x;} 

22 : return T 


A node can be added several times to Narrow during FIM execution, since 
every time an upwind (parent) neighbor is updated, the node can improve its 
value. In the worst case. Narrow contains the whole grid and the loop would go 
through all the points several times. Operations on the list are 0(1). Therefore, 
the overall computational complexity of FIM is 0{n). 

For a small enough e (depending on the environment), FIM will return 
the same solution as FMM. However, it can be speed up allowing small errors 
bounded by e. 
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6 Experimental Comparison 

6.1 Experimental setup 

In order to give an impartial and meaningful comparison, all the algorithms 
have been implemented from scratch, in C++11 using the Boost.Heap library"^. 
An automatic benchmarking application has been also created so that the 
experiments are carried out and evaluated in the most systematic possible 
way. 

This implementation is focused on time performance and compiled using 
G++ 4.9.2 with optimizations flag -Ofast. However, no special optimizations 
have been included. All algorithms use the same primitive functions for grid 
and cell computations. The times reported correspond to an Ubuntu 14.04 64 
bits computer running in a Dual Core 3.3 GHz with 4Gb of RAM. However, all 
experiments were carried out in one core. Only propagation times are taking 
into account. The computation time used in the initialization has been omitted 
since it can be done offline, besides it is similar for all algorithms and represents 
a little percentage of the total computation time. 

Although the algorithms are deterministic (their result only depends on 
the grid, which is not modified), the results shown are the mean of 10 runs for 
every algorithm, so that the deviation of the results is practically 0. 

For UFMM, the default is a maximum range of T of 2 units and 1000 
buckets (the checkerboard experiment required different parameters, see Sec¬ 
tion 6.1.4). The e parameter for FIM is set to 0 (actually 10“^^ to provide 
robust 64bit double comparison). 

Although error analysis is not in the scope of this paper, these can be com¬ 
pared among the existing papers since they are implementation-independent. 
UFMM errors are reported in those experiments with non-constant velocity. 
Usually, Li and L^o norms of the error are reported. Most of the works com¬ 
pute norm Li as: 


irii = 5 ] iTii (18) 

xiea” 

where X is treated as a regular vector. However, following [7], we treat the 
numerical solutions as elements of Lp spaces (generalization of the p-norm to 
vector spaces), where Li norm is defined as an integral over the function. The 
result is a norm closely related to its physical meaning and independent of the 
cell size. Li is numerically integrated over the domain and therefore computed 
as (assuming cubic cells and grids): 

irii = mh^\ = i^ii (19) 

Four different experiments have been carried out, which represent charac¬ 
teristic cases. By combining these problems it is possible to get close to any 


The source code is available at https://github.com/jvgomez/fastmarching 



Fast Methods for Eikonal Equations: an Experimental Survey 


21 


situation. These experiments were chosen attending also to the most common 
situations tested in the literature. 


6.1.1 Empty map 

This experiment is designed to show the performance of the methods in the 
most basic situation, were most of the algorithms perform best. An empty 
map with constant velocity represents the simplest possible case for the Fast 
Methods. In fact, analytical methods could be implemented by computing 
the euclidean distance from every point to the initial point. However, it is 
interesting because it shows the performance of the algorithms on open spaces 
which, in a real application, can be part of large environments. 

The same environment is divided into a different number of cells to study 
how the algorithms behave as the number of cell increases. Composed by an 
empty 2D, 3D and 4D hyper-cubical environment of size [0,1]'^, with N = 
2, 3,4. Constant velocity F; = 1 on f2. The wavefront starts at the center of 
the grid. 

The number of cells was chosen so that an experiment has the same (or as 
close as possible) number of cells in all dimensions. For instance, a 50x50 2D 
grid has 2500 cells. Therefore, the equivalent 3D grid is 14x14x14 (2744) and in 
4D is 7x7x7x7 (2401). This way, it is possible to also analyze the performance 
of the algorithms for a different number of dimensions. Thus, we have chosen 
the following number of cells for each dimension for 2D grid: 

2D : {50,100,200,400,800,1000,1500,2000,2500,3000,4000} 

Consequently, the 3D and 4D cells are: 


3D :{14, 22, 34, 54, 86,100,131,159,184, 208, 252} 
AD :{7,10,14,20,28, 32, 39,45,450, 55,63} 


6.1.2 Alternating barriers 

In this case, we want to analyze how the algorithms behave with obstacles 
(Fi = 0) in a constant velocity environment (Fi = 1). The obstacles cause the 
characteristics to change. 

The experiment contains a 2D environment of constant size [0, l]a;[0, 2] 
discretized in a 1000x2000 grid. A variable number of alternating barriers 
are equally distributed along the longest dimension. The number of barri¬ 
ers goes from 0 to 9. Examples are shown in Fig. 4. Analogously, in 3D a 
[0, l]a:[0, l]a;[0, 2] environment represented by a 100x100x200 is chosen, with 
equally-distributed alternating barriers (from 0 to 9) along the z axis. The 
wavefront starts in all cases close to a corner of the map. 
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(a) 1 barrier. (b) 5 barriers. (c) 9 barriers. 

Fig. 4 2D alternating barriers environments. 


6.1.3 Random velocities 

In order to test the performance of the algorithms with random velocities, or 
noisy images as in the case of medical computer vision. This experiment creates 
a 2D, 3D and 4D environment of size [0,1]^ with N = 2,3,4 discretized in 
a 2000x2000 grid in 2D, 159x159x159 in 3D and 45x45x45x45 in 4D. These 
discretizations are chosen so that it is possible to compare directly with the 
empty map problem for the corresponding grid sizes. The wavefront starts in 
the center of the grid. 

Additionally, the maximum velocity is increased from 10 to 100 (in steps 
of 10 units) to analyze how the algorithms behave with increasing velocity 
changes. 2D examples are shown in Fig. 5. 



(a) Max. velocity = 30 (b) Max. velocity = 60. (c) Max. velocity = 100. 

Fig. 5 2D random velocities environments. Lighter color means faster wave propagation. 


6.1.4 Checkerboard 

The random velocities experiment already tested changes in velocities. How¬ 
ever, those are high-frequency changes because it is unlikely to have 2 adjacent 
cells with the same velocities. In this experiment low-frequency changes are 
studied. The same environment and discretization as in random velocities is 
now divided like a checkerboard, alternating minimum and maximum veloci¬ 
ties. Analogously, the maximum velocity is increased from 10 to 100, while the 
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minimum velocity is always 1. There are 10 checkerboard divisions on each 
dimension. The wavefront starts in the center of the grid. 2D examples are 
shown in Fig. 5. 



(a) Max. velocity = 30 (b) Max. velocity = 60. (c) Max. velocity = 100. 

Fig. 6 2D checkerboard environments. Lighter color means faster wave propagation. 


In this case, UFMM in 3D and 4D performed very poorly with the de¬ 
fault parameters. Our previous experiments show that the best parameters for 
UFMM are approximately 1000 buckets with a maximum range of 0.01 in 3D. 
In the 4D case, 20000 buckets and maximum range of 0.025 are used. 


6.2 Results 
6.2.1 Empty map 

An example of the times-of-arrival field computed by FMM is shown in Fig. 7. 
Note that all algorithms provide the same exact solution in this case. The 
higher the resolution the better the accuracy. 



(a) 50x50 (b) 800x800 (c) 4000x4000 


Fig. 7 Example of the resulting times-of-arrival maps applying FMM to the empty envi¬ 
ronment in 2D. 
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The results for the empty map experiment are shown in Fig. 8 for 2D, Fig. 9 
for 3D, and Fig. 10 for 4D. In all cases 2 plots are included: raw computation 
times for each algorithm, and time ratios computed as: 


Alg. Time 
ratio = ^ ^ r^ — 

FMM Time 


( 20 ) 


In all cases, both LSM and DDQM are the fastest algorithms. DDQM tends 
to perform better in smaller grids, specially en 4D. This was an expected result 
since this is the case in which they perform less sweeps (queue recomputations 
in DDQM). Besides, FSM is slower than UFMM in all cases and than FIM in 
3D and 4D. As velocities are constant, UFMM provides the same solution as 
other methods, 

GMM slightly improves FMM and FMMFib in 2D. However in higher di¬ 
mensions this improvement becomes larger, but it always performs worse than 
UFMM, FIM, DDQM and LSM. In a previous comparison between GMM and 
FMM [19], GMM was about 50% faster in all cases than FMM. In this results 
GMM is at most 40% better. We attribute this difference to implementation, as 
the heaps for FMM and FMMFib are highly optimized. FIM is always faster 
as it only needs one iteration trough the narrow band, while GMM always 
performs two. 

SFMM results are of special interest since it is a minor modification of 
FMM but, however, highly outperforms FMMin all cases, and even FIM in 2D 
and small 3D grids. In 4D FIM becomes faster. Besides UFMM outperforms 
FSM. 

As expected, FMMFib is worse than FMM for almost all sizes. However, 
when dimensions increase FMMFib quickly outperforms FMM as the number 
of elements in the narrow band increases exponentially with the number of 
dimensions, therefore the better amortized times of Fibonacci heap become 
useful. 



Fig. 8 Computation times and ratios for the empty map environment in 2D. 
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# Cells ^,(,7 # Cells 


(a) Computation times. (b) Time ratios against EMM. 

Fig. 9 Computation times and ratios for the empty map experiment in 3D. 



# Cells 

(a) Computation times. 


(b) Time ratios against FMM. 


Fig. 10 Computation times and ratios for the empty map experiment in 4D. 


6.2.2 Alternating harriers 

An example of FMM results in some of the alternating barriers environment 
is shown in Fig. 11, while performance results (times and ratios) for 2D and 
3D are shown in Fig. 12 and Fig. 13 correspondingly. 

In this case, results are very similar to the empty map experiment. However, 
as the number of barrier increases and the environment becomes more complex, 
FSM and LSM decrease their performance. However, LSM is still faster than 
some algorithms in most cases. 

DDQM also suffers from map complexity, however it affects much less and 
only in 3D. Changes in characteristic directions require more sweeps to com¬ 
pute the final times-of-arrival map. Therefore, FSM and LSM become worse 
as the number of barriers increases, despite the fact that the more barriers the 
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less cells are to evaluate. This is the reason why all other algorithms tend to 
lower times as the number of barriers increases. 

UFMM provides again the same solution as other methods. It is the second 
fastest algorithm behind DDQM, closely followed by SFMM and FIM. 



(a) 1 barrier. (b) 5 barriers. (c) 9 barriers. 


Fig. 11 Example of the resulting times-of-arrival maps applying EMM to some of the 
alternating barriers environment in 2D. 



Fig. 12 Computation times and ratios for the alternating barriers experiment in 2D. 


6.2.3 Random velocities 

The output of the FMM for the random velocities map is apparently close to 
the one of the empty map (Fig. 14), but with the wavefronts slightly distorted 
because of the velocity changes. 

Although the results on the times-of-arrival map are slightly different from 
the results obtained in the empty map, the performance of the algorithms 
are highly modified. 2D, 3D and 4D results are shown respectively in Fig. 15, 
Fig. 16, and Fig. 17. In this case, raw computation times are shown together 
with a zoomed view of the fastest algorithms to make the analysis easier. Note 
that all methods become slower with non-constant velocities. The reason is 
that the narrow band tends to have more elements in this cases. 
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4 6 

# Barriers 


(a) Computation times. 


(b) Time ratios against FMM. 


Fig. 13 Computation times and ratios for the alternating barriers experiment in 3D. 


Some algorithms become unstable with non-uniform velocities: FSM, LSM 
and DDQM. DDQM is able to maintain the fastest time for slight velocity 
changes. But when these are sharper the double-queue threshold becomes un¬ 
stable. However, for a high number of dimensions this effect vanishes (its com¬ 
putation time is barely modified with the number of dimensions) and DDQM 
becomes the fastest algorithms together with FIM and GMM. 

SFMM provides one of the best performances across all dimensions. And 
FIM becomes relatively faster in 3D and 4D. However, in most of the cases 
GMM is the fastest algorithm. FIM requires now multiple iterations to con¬ 
verge to the solution while GMM guarantees convergence with only 2 itera¬ 
tions. 

Finally, UFMM requires special attention as it does not return the same 
solution than the other Fast Methods. Its performance is highly affected by 
the number of dimensions. The main reason is the election of the parameters: 
they were experimentally chosen to optimize 2D performance. However, these 
parameters are no longer useful in other number of dimensions. We consider 
UFMM parameter tuning to be a complex task. 
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(a) Max. velocity=30 (b) Max. velocity=60 (c) Max. velocity=100 


Fig. 14 Example of the resulting times-of-arrival maps applying FMM to the random 
velocities environment in 2D. 
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Fig. 15 Computation times and ratios for the random velocities experiment in 2D. 




Fig. 16 Computation times and ratios for the random velocities experiment in 3D. 


Table 2 summarizes the largest errors for this experiment. As the number 
of dimension is increased, the error decreases exponentially while the compu¬ 
tation time increases exponentially. Therefore, by properly tuning parameters 
for 3D and 4D better times could be achieved while keeping a negligible error 
in most practical cases. 


Table 2 Largest Li and Loo errors for UFMM in the random velocities experiment. 
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(a) Computation times. 



40 60 

Max. Velocity (units) 


(b) Zoomed view. 


Fig. 17 Computation times and ratios for the random velocities experiment in 4D. 


6 . 2.4 Checkerboard 

Finally, different times-of-arrival map returned by FMM applied to the checker¬ 
board map are shown in Fig. 18. Numerical results of the computation times 
are included in Fig. 19 for 2D, Fig. 20 for 3D, and Fig. 21 for 4D. 

The results are relatively close to the random velocities experiment. How¬ 
ever, the differences for FSM, LSM and DDQM are much smaller. In fact, 
DDQM presents a poor performance in 2D, but in 3D and 4D it becomes the 
fastest algorithm for higher velocity modifications. 

GMM and FIM have very similar results for 3D and 4D. Since the envi¬ 
ronment is structured, FIM does not require too many iterations to compute 
the final value. 



(a) Max. velocity=30 (b) Max. velocity=60 (c) Max. velocity=100 

Fig. 18 Example of the resulting times-of-arrival maps applying FMM to the checkerboard 
environment in 2D. 


UFMM errors are shown in Table 3. In this case, UFMM becomes worse 
with the number of dimensions but it is among the fastest algorithms in all 


cases. 
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Fig. 19 Computation times for the checkerboard experiment in 2D. 




Fig. 20 Computation times for the checkerboard experiment in 2D. 


The differences between this experiments and random velocities is that the 
environment presents a well-defined structure, and locally acts as a constant 
velocity environment, as in empty map experiment. 


Table 3 Largest L\ and L^o errors for UFMM in the checkerboard experiment. 
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7 Discussion 

With the four experiments designed, the main characteristics of the Fast Meth¬ 
ods have been shown. Any other environment can be thought as a combina- 
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(a) Computation times. 


(b) Time ratios against EMM. 


Fig. 21 Computation times for the checkerboard experiment in 4D. 


tion of free space with obstacles and high-frequency or low-frequency velocity 
changes of different magnitude. Note also that the grid sizes covered by the 
experiments vary from extremely small to extremely big grids. In practice, it 
is hard to find applications requiring more than 16 million cells. 

FIM results can be speeded up for non-constant velocity problems if larger 
errors are allowed. UFMM can be probably improved as well. However, our 
experience is that the configuration of its parameters is complex and requires 
a deep knowledge of the environment to be applied on. 

Several conclusions can be extracted from the conducted experiments: 1) 
There is no practical reasons to use FMM or FMMFib as SFMM is faster 
in all cases with the same behaviour as its counterparts. 2) If a sweep-based 
method is required, LSM should be always chosen, as it greatly outperforms 
FSM. 3) In problems with constant velocity DDQM should be chosen, as it has 
shown the best performance for the empty map and the alternating barriers 
environments. 4) For variable velocities, but simple scenarios, GMM is the 
algorithm to choose. 5) UFMM is hard to tune and the result has errors. 
Also, it has been outperformed in most of the cases by DDQM in constant 
velocity scenarios, or by SFMM or FIM in experiments with variable velocities. 
6) There is not a clear winner for complex scenarios with variable velocity. 
UFMM can perform well in all cases if tuned properly. Otherwise, SFMM is a 
safe choice, specially in cases where there is not too much information about 
the environment. 

If a goal point is selected, cost-to-go heuristics could be applied [36], and 
thus enormously affect the results. Heuristics for FMM, FMMFib and SFMM 
are straightforward. They would improve the results in most of the cases. They 
could also be applied to UFMM. However, it is not clear if they can be applied 
to other Fast Methods. It is also of interest the solution of anisotropic problems 
[8], solved only by FMM-based methods. 
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8 Conclusions 

Along this paper we have introduced the main Fast Methods in a common 
mathematical framework, focusing on a practical point of view. 

With this work we aimed at closing the discussion about which method 
should be used in which case, as this is the first exhaustive comparison of the 
main Fast Methods (up to our knowledge). 

The code is pnblicly available as well as the automatic benchmark pro¬ 
grams. This code has been deeply tested and it can serve as a base for futnre 
algorithm design, as it provides all the tools required to easily implement and 
compare novel Fast Methods. 

The future work focuses in 3 different aspects: develop the analogous 
work for parallel Fast Methods [9], study the application of these methods 
to anisotropic problems and also to the new Fast Marching-based solutions 
focused on path planning applications [18,38]. Finally, the combination of 
UFMM and SFMM seems straightforward and it would presumably outper¬ 
form both algorithms. 


Acknowledgements Authors want to gratefully acknowledge the contribution of Pablo 
Gely Munoz to the GMM, FIM and UFMM implementations, and to Adam Chacon for the 
interesting discussions and suggestions towards improving the work. 
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